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A family of models of liquid on a 2D lattice (2D lattice liquid models) have been proposed as 
primitive models of soft-material membrane. As a first step, we have formulated them as single- 
component, single-layered, classical particle systems on a two-dimensional surface with no explicit 
viscosity. Among the family of the models, we have shown and constructed two stochastic models, 
a vicious walk model and a flow model, on an isotropic regular lattice and on the rectangular 
honeycomb lattice of various sizes. In both cases, the dynamics is governed by the nature of the 
frustration of the particle movements. By simulations, we have found the approximate functional 
form of the frustration probability, and peculiar anomalous diffusions in their time-averaged mean 
square displacements in the flow model. The relations to other existing statistical models and 
possible extensions of the models are also discussed. 

PACS numbers: 87.16.dj, 87.10.Hk, 87.16.aj, 05.50.-fq 
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I. INTRODUCTION 

Two-dimensionally extended materials, or membranes, 
have been paid much attention to in various fields of 
study, such as lipid bilayers and surfactants Q in 

soft- matter and biophysics . Alongside, many theoret- 
ical studies on those membrane-type objects have been 
done toward a general understanding of physical proper- 
ties of such materials 0, @, 0| • Some of them were simu- 
lation approaches iS,j6| and others were with field theo- 
retical techniques [l|, Q • For example, one of the authors 
studied the pore formation in a binary giant vesicle of 
two-component lipid bilayer, and reproduced its dynam- 
ical shapes in the field-theoretical approach with some 
numerical calculations [1]. Although most of the works 
contributed successfully to our knowledge of individual 
membranes, it can be said that we are yet very far from 
the general understanding. 

In this work, we aim to give some new insight to the 
general understanding of the physical properties of soft- 
material membranes by studying mathematical models 
on discretized surfaces both analytically and numerically. 
That is to say, we propose a family of fiuid particle mod- 
els of two-dimensional lattice as simpler models of mem- 
branes: tightly pack molecular particles on a given lat- 
tice, and give the particles' dynamics by some stochastic 
process, or by a set of mathematical rules. For more 
simplicity, we consider single-component, single-layered 
incompressible systems, and assume no explicit viscosity 
and an upper bound for the fiow velocities on the parti- 
cles. After the formulation, we construct two stochastic 
models of 2D lattice liquid, as our first steps. One of 
them is found to be the vicious random walk Q in its 
densest phase and a version of the dynamic lattice liquid 
model (DLL model) [9j] in two dimensions. The dynam- 
ics of these models is given by completely independent 
time steps and would not generate various fluid dynam- 
ical phenomena such as lasting convective flows, while 
our formulation does not exclude such possibilities. In 
fact, the other model is constructed as the fiow model 



with a mathematical rule, which has no corresponding 
known models. We further restrict ourselves to a finite 
flat honey-comb lattice and simulate the two models with 
their highest mobility. Comparing with similar lattice 
models such as various lattice random walks or the loop 
gas model (SOS model) we show some statistical 

properties of the models by the simulation results. 

The results tell that the flow model may exhibit some 
anomalous diffusion, in terms of time-averaged mean 
square displacements (TAMSD), and the ergodicity with 
its fairly long relaxation time. Relations to other models 
and some possible extensions of the models will also be 
discussed in the final section. 



II. PRELIMINARIES 

There are various situations of the membrane systems, 
and each situation would deserve a single mathemati- 
cal model in order to simulate it. However, since we 
aim to construct a basic model to be utilised for those 
various situations, we focus exclusively on a very sim- 
ple but basic situation among them as a building block. 
First, we focus on a single-component and single-layered 
molecular system on a two-dimensional surface, treating 
the surrounding medium as a kind of heat bath to the 
surface. In other words, we introduce some rules or ex- 
ternal forces onto the elements of the membrane rather 
than consider the interactions between the membrane 
and the medium. Secondly, we focus exclusively on a 
simple classical interaction such as "collision" , neglecting 
any other complicated interaction between the elements 
of the membrane. Let us call this 'no explicit viscosity'. 
Thirdly, we assume that the density fluctuation can be 
negligible and two elements would not occupy the same 
location of the surface at any time. This would corre- 
spond to the incompressibility condition of the surface 
from the two-dimensional fluid point of view. To sum- 
marise, we assume and focus on a single-component and 
single-layered incompressible particle systems with no ex- 



2 



plicit viscosity on a two-dimensional surface exposed in 
an external (random) field. 

Now, we are ready to define the 2D lattice liquid mod- 
els. Given a two-dimensional lattice A, a 2D lattice liq- 
uid model can be defined as a mathematical model of 
the particles which can only reside at the sites of the lat- 
tice A, covering all or almost all of the sites, and which 
flow inside A as molecular-dynamical units of liquid do 
inside them. Here, the particles represent the elements 
of the membrane, and are assumed to be of the same 
physical properties: their sizes, masses, and so on. The 
word 'liquid' comes directly from two parts of the above 
definition, one of which is that the particles cover all or 
almost all of the sites of A, and the other of which is that 
the particles flow. Obviously, there is an ambiguity on 
the phrase "all or almost all". To avoid this ambiguity 
in this paper, we force each site to hold only and exactly 
one particle at a time, covering all the sites. We call it 
"the occupation number is exactly one everywhere" . Al- 
though we can make some of the sites empty or put as 
many particles as possible, the assumption of the density 
fluctuation eliminates the latter and we simply omit the 
former in this paper. Note that the definition of the lat- 
tice liquid would not exclude these cases and would be 
worth studying elsewhere. 

One can consider the boundary conditions in the pres- 
ence of boundary of the surface. There are two types of 
boundary conditions, one of which is on particles' dynam- 
ical degrees of freedom, and the other of which is on their 
internal degrees of freedom, such as particle species. The 
former type of boundary conditions are naturally embed- 
ded in the lattice on which we define the model. On the 
latter, since we only consider single-component systems, 
there is no boundary condition on the particle species. 

To complete the definition of the model, the particles' 
"flows" should be defined by imposing a set of mathe- 
matical rules on the flows. This consequently determines 
the dynamics of the system. Before going into the details 
of the rules, let us introduce the lattice conventions and 
consider the configuration space of the model to begin. 



A. The configuration space of tiie model 

A lattice A is defined by a set of sites i G V{A) and a 
set of bonds ej j- G E(A). A bond e^j is connecting a pair 
of the sites, i,j € ^(A), which are called neighbouring 
sites or in the neighbourhood. The lattice spacing a - the 
shortest length of the bonds specifies a fine aspect of the 
lattice. We consider only the cases with the same bond 
length, unless otherwise stated. For that reason we set 
a = 1 without loss of generality. The total number of the 
sites and the total number of the bonds are denoted by 
dim{V{A)) and dim{E{A)), respectively. We abbreviate 
them by dimV and dimE for simplicity. The degree of 
the site, deg{i), is the number of the bonds connected to 
the site i. If there is only one value for deg{i) on A, we 
write it by deg for short, and we consider mainly such 



cases. In the later sections of some numerical calcula- 
tions, we focus on the honeycomb lattices of various sizes 
as a series of examples. 

On a given lattice A, under any set of the rules of the 
movements, a flow of a particle can be sketched schemat- 
ically in FigU] below. Introducing the minimum time 
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FIG. 1: One of the "world-lines" illustrating a particle trav- 
elling on a lattice. 

lapse (5t > 0, a particle at the site i at the time t flew 
from the site (i — Vi) at t — St with the velocity Vi, where 
St is implicitly taken as a unit {St — I). Then, the par- 
ticle at i at t flows toward {i + Vi + Svi) at the future 
time t -\- St, subject to some fluctuation Svi by exter- 
nal forces. Explicitly writing the time dependence of the 
variables, Vi{t) denotes the flow velocity of the particle 
from t — St to t, and Vi{t) + Svi{t) does from t to t + St. 
We denote the latter velocity by a redundant variable 
v'j^{t) = Vi{t) + Svi{t) for convenience. The solid con- 
tinuous line in FiglT] represents the particle travelling in 
(2-|-l)-dimensional space-time, and all the particle move- 
ments can be represented by such world-lines. Namely, 
the conflguration of the system at a certain time t, can 
be uniquely given by a set of the flow velocities: 

{v^{t)) for all i e V{A). (1) 

Accordingly, the time evolution is realised by mapping a 
configuration {vi{t)} to another configuration {vi{t+St)} 
at a successive time, while a particle at i travels to i' with 
the velocity v'^{t) such that: 

i'^i + v'S), v,,{t + St)=v[{t). (2) 

For consistency, Vi and v[ cannot take completely arbi- 
trary values, but are subject to the constraint that {vi\ 
and {v'j} map V{A) to itself and they are bijective. In 
other words, the configuration space is given by the whole 
set of such bijections. Such a space is highly complicated 
in general, however one assumption actually makes the 
space more tractable. 

We assume the following upper bound for the veloci- 
ties: 

\v,\ < 1 for afl i e V{A). (3) 
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This means that the particles are prevented from jump- 
ing in between non-nearest neighbouring sites. There- 
fore, the assumption is not too restrictive in the physics 
context. When the bond lengths are all the same, this 
leads the absolute value \vi\ takes either or 1. Accord- 
ingly, the velocity Vi takes one of {deg{i) + 1) values. 
Since \v'^\ obeys the same, the velocity change 5vi takes 
at most one of {deg{i) + 1) values as well. 

With the upper bound of the velocities and the as- 
sumption of the occupation number, a configuration of 
the particles can be represented by a collection of ori- 
ented non-intersecting loops on A. In order to illustrate 
this fact, let us give a set of random velocities, Svi{t), 
with the condition Vi{t) = 0. All particles potentially try 
to move toward the directions they have been kicked by 
some external forces. However, this does not necessarily 
give movements to the particles due to the frustration of 
such momenta (Figl^Ja)). The frustration is an outcome 
of the excluded volume effect of the particles together 
with the restriction of the occupation number. In order 
to generate convective flows, some part of the "potential" 
movements generated by the forces must be enclosed so 
that the particles can move (Figl^fb)). Hence, the rel- 




(a) (b) 

FIG. 2: (a) An example of the frustration of the velocities 
on a triangular lattice. Both red and blue sequences of vec- 
tors show the frustration, (b) An example of the possible 
movements on a triangular lattice. Red and blue sequences 
of vectors show the circulations. 

evant part of the configuration {vi} can be represented 
by a collection of oriented non-intersecting loop graphs 
on A. We label the space of the oriented non-intersecting 
loop configurations by C{A). We also label the loop con- 
figuration without the orientations by C{A) for later use. 
Formal definition of the (oriented) loop configurations in 
graph theory, refer to Appendix [Xj Note that, if and only 
if Vi{0) = 0, the above algorithm coincides with a ver- 
sion of the DLL model In what follows, we focus on 
the situations where the loop configuration solely plays 
a crucial role. 



B. Formulation of the model 

As was stated before, the set of mathematical rules of 
the particle-flows give the principal part of the definition 
of the model. However, regardless of how and what sort 
of the rules are imposed, they are to be expressed in the 
form of the velocities and the velocity changes: {vi} and 
{5vi}, or in the form of the present loop configuration 



and the resulting one, G, G' G C{A), representing the ve- 
locities. So, we formulate the 2D lattice liquid models as 
a family of statistical models defined with such variables. 

One of the most straightforward expressions of the 
propagator would take the form of: 

g{xi,x2,--- ;G,t;G',t') 

or g{xi,X2,--- ;{v^it)},t■,{v,{t')},t'), (4) 
or as its ensemble: 

ZS{xi,X2,--- ;G,t;t') 

= E G{xi,X2,--- ■,G,t;G',t'), (5) 

G'e£(A) 

where the x^'s are other unspecified parameters to pre- 
scribe the propagator. If all the processes are statistical 
from the beginning, i.e., even the initial configuration is 
given by some distribution, the following quantity would 
be the partition function in question: 

= J2 Prmt{G') ZS{xi,x2,--- ;G,t;t'), (6) 

Ge£(A) 

where the function Pmit (S) is the statistical weight of the 
configuration S. Practically, we will focus on the second 
expression ([5]) of the propagator and the first one (|4|) for 
the construction of the second. 

Let us construct a more concrete yet abstract form of 
the propagator by the probabilities of the configurations. 
For the time step of St, given the present configuration 
G, let us define p{i) as the "potential" probability of a 
particle remaining at the site i, and x{i;j) is the "poten- 
tial" probability of the particle at i being kicked toward 
one of the nearest neighbouring sites j. They satisfy: 

p{z)+J2x{i;j)^l, (7) 

j 

where the summation is over all nearest neighbouring 
sites to i. In general, p{i) and x{i;j) may carry the infor- 
mation from the past, but we ignore it in the rest of this 
manuscript. Recall that some of the kicks with the "po- 
tential" probabilities might not be performed due to the 
before mentioned frustration. Suppose non-frustrated 
part of the particles are expressed in the configuration 
G' . The complementary graph can be given with its pro- 
jection G' on C{A) by G" = G'. Then, the probability of 
having G' in A, P{G'), can be given by the product of the 
loop probability Pi of having G' in G' and the probability 
Pf of the frustration in G' as follows: 

P(G')-P,(G')-P/(G'), (8) 

where the initial configuration G is implicit, and the loop 
probability can be expressed by multiplying the proba- 
bilities x{i;jys along the oriented loops: 

PiiG')^ n ^(*;^')- (9) 
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Note that the definitions of the above probabilities need 
some minor extensions to the subgraphs for their con- 
sistency. For the detailed definition of the probability 
fimctions and their properties, refer to Appendix iBl 

The probability P{G') is nothing but the first expres- 
sion (ID) of the propagator for one time step. Accordingly, 
the second expression ([SJ of the propagator for one time 
step is: 

ZS{G,t;t + 6t) = J2 ^(^')- (10) 

G'G£(A) 

Specifying its G dependence by P{G'; G), the propagator 
for a longer time span can be expressed by: 

Z9{G^°\t;t') 

= E ■•• E flP{&'^;G^'-'^)An) 

G(i)6£(A) G(-ff-i)e£(A) 

where t' = t+K6t. Hence, the construction of the propa- 
gator reduces to a detailed construction of the probability 
functions, Pj and Pi. Note that the above expression ac- 
tually corresponds to the discretised version of the path 
integral. Given the probabilities, an expectation value of 
some physical quantity O for St can be given as usual by 
the form of: 

(0)c5 = E 0-P{G';G). (12) 

G"e£(A) 

III. MODELS 

It is obvious that the set of the mathematical rules 
of the particle-flows give the principal part of the def- 
inition of the 2D lattice liquid model as well as of its 
formulation. In fact, one can easily imagine the follow- 
ing three sorts of the rules for the particle-flows. Namely, 
the particles flow according to some stochastic process, or 
subject to some external thermodynamic forces from the 
surrounding medium, or obeying some physics laws with 
their (fluid) dynamical properties such as the Langevin- 
type equation or the Navier-Stokes equation in the low 
T single-component limit. In order to avoid ambiguities 
on discretisation and limiting procedures in the equation- 
based approaches, we deal only with the stochastic cases 
here. Some relations to thermodynamic or equation- 
based approaches will be discussed later. 

A. Stochastic Markovian model (I) 

In this section, we construct and formulate one of 
the simplest stochastic models, the stochastic Markovian 
model, whose rules are as simple as assigning simple kicks 
to push the particles and arresting them at the end of 
each step. The model is turned out to be equivalent to 



the vicious random walk in its densest phase and a ver- 
sion of the DLL model on two-dimensional lattice. The 
precise definition of the rules are given as follows. 

First, let p{i) be non-negative constant, and x{i]j) be 
of isotropic kicks such that: 

p{i) ~ p for < p < 1, 

Here, we restrict them to be strictly independent of the 
past, ensuring that the process is Markovian. Second, 
we assign the kicks according to the above probabilities, 
p and x(i), and those kicks are projected onto some ori- 
ented loop configuration G. Perform the moves of G. 
Thirdly, as a key part of the rules, we force: 

Vj{t + dt) ^6v,{t) for alH, (14) 

where j = i + Svi{t). This means that the configuration 
of the velocities is set at every step by the fiuctuation 
{Svi} only. Accordingly, one of the two sets, {vi} and 
{Svi}, becomes redundant. Note that \Svi\ inherits the 
upper bound from \vi\: \6vi\ < 1. Practically, it follows 
1^-^,1=0,1. 

The probability of an oriented loop configuration G on 
the sublattice G can be written by: 

PiiG)= n ^(*)- (15) 

e,.je£;(G) 

If (ieg(i)'s are all the same as of regular lattices, the above 
probability can be simplified as: 

Pi{G;x)=x'^'^\ (16) 

where we specify x as an argument, x = {l—p)/{deg), and 
1{G) denotes dimE{G) = \6vi\ for {6vi} of G, i.e., the 
total length of the loops in G. Note that x plays the role 
of the fugacity of the bonds in G. If we further define 
the probability of a non-oriented loop configuration G on 
the sublattice G, summing up all the orientations of G, 
it follows that: 

P,(G;x) = 2"(<=)a;'(^\ (17) 

where Z(G) — 1{G) by definition, a{G) counts the num- 
ber of the loops in G, and the factor '2' indicates the 
directions of each oriented loop in G. 

Accordingly, the propagator, or the partition function, 
can be written down by: 

g''°{x)= Pi{G;x)Pf{G;x), (18) 

Gez:(A) 

which is equivalent to unity by definition. Pf{h.\ x) gives 
the probability of the no-loop configuration. The prob- 
abilities Pf and Pi take physical values, < 1 and 
|Pi I < 1 for < p < 1. It should be mentioned here that 
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if Pf (G) is set to be one for any G then the partition func- 
tion reduces to that of the non-oriented non-intersecting 
loop configurations on a lattice, which is known as the 
loop gas model [l^: Z'°°P(n,a;) with n = 2. So, P/(G) 
plays a crucial role in our model. 

As in the definition (|14|) . the present configuration is 
always trivial so that the above expression (IT51) is exact 
at any t. Therefore, the propagator can be expressed by: 



(19) 



Hence, it can be summarised and simplified to the prop- 
agator of the time lapse t as follows: 



t/5t 



(20) 



The expectation value of a physical quantity O is given 
in a conventional manner: 



{0)t 
(O) 



J2 0-Pf{G;x)Pi{G- 

Gec{A) 



(21) 



Before giving our simulation results, let us show some 
simple relations and calculations on an isotropic regular 
lattice. 



1. On an isotropic regular lattice 

On an isotropic regular lattice, such as a square lattice 
with periodicities, the one-point function of the velocity 
can trivially be given by: 



0, 



(22) 



by symmetry. Similarly, the expectation value of the ab- 
solute value of the velocity at i can be given by: 



dimV 



I 



dimV 



for all i, 



(23) 



where Tis the expectation value of the length of the loops: 
I = X]gg£(A) KG)Pf{^i x)Pi{G; x). Besides, the variance 
of the velocity becomes: 



xAv^r) = {\v,\^) = 



I 



dimV 



(24) 



where the normalised lattice spacing is used. 

With the above facts, one can state that each particle 
behaves like a random walker on the given lattice with 
the constant trapping ratio of p = 1 — ( where a 
particle moves to a certain direction with the probability 
X = ( \ vi\)/deg. Namely, for any to for large A or small 
t: 



{xj{to+t)-xi{to)) = 0, 



[ixi{to+t)-Xi{to)f) 



X ■ deg 



(25) 



where xi{t) denotes the position of the /-th particle at 
the time t. However, it should be noted that the system 
is not a collection of independent simple random walks, 
as our particles do not occupy the same site at the same 
time. In fact, the system is in the densest phase of the 
vicious random walks — mutually avoiding N particles' 
random walks on a lattice with N sites Q . The difference 
between this system and a collection of the independent 
simple random walks can be clearly shown by the spacial 
correlation of the velocities. For example, the spacial 
correlation of a nearest pair of the system is given by: 



\Vi\\V, 



> 



21 



dimV ■ deg 



for 



(26) 



whereas those of the simple random walks is given by a 
simple product of one-point functions: < >=< 

1^,1 X \v,\ {^y. 



B. Stochastic flow model (II) 

In the previous sections, we have formulated the model 
where the streamlines are temporally generated and live 
only for the duration of 6t. Because the probability of 
the shortest loop in a honeycomb lattice is at most 
in the bulk, I seems to be very small. We infer in the 
model (I) that the stochastic process realises the solid- 
like nature of the membrane rather than its liquid nature 
as the DLL model targeted supercooled liquids. In this 
section, we intend to explore another class of stochastic 
models where the particles keep flowing on the membrane 
rather than obeying the 'stop and go' evolution. In order 
to realise this phenomenon, we take the simplest way: 
if no fluctuations are added to the system, the particles 
keep flowing along the oriented loops. This rule can be 
expressed by the relation in the absence of the fluctua- 
tions: 



v,{t + St) = Vi{t). 



(27) 



Now, it is evident that the key point of the construction 
lies on how to add stochastically generated fluctuations 
onto the existing streamlines. 

In order to avoid complications for realising the rela- 
tion (P7)) and the frustrations in terms of {5vi}, let us 
formulate the above addition in terms of the oriented 
loop configurations. Say, Gi is the initial configuration 
and G/ is the resulting configuration. The arithmetic 
rule (addition/subtraction) for the oriented loop config- 
urations can be illustrated in a multiplicative form by: 



^G • G, = G/ e C{A), 



(28) 



where SG stands for the stochastically generated fluctu- 
ations, and vanishing SG is given by the no-loop conflg- 
uration /(A) naturally satisfying the rule ([?7|) . In C{A), 
there is a natural definition of such an additive operation 
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as the symmetric difference. We extend and modify it to 
our case for C{A). 

Let Gi, G2 be eiements in C{A). If G2 consists of the 
same sites and bonds with opposite directions as those 
of Gi, then we call G2 the inverse of Gi, denoting it by 
G2 ~ Gi^. The additive operation between Gi and G2 
is defined in the multiplicative form by: 

Gi • G2 = G, (29) 

such that the directed bonds of the elements satisfy: 

G = GiUG2-GinG2"^-Gj;inG2, (30) 

while the disconnected sites are removed from G. The 
above addition trivially satisfies the followings: 

Gi • /(A) = /(A) • Gi Gi for aU Gi e £(A), 
Gi-G2 - /(A) only when G2 = G^S (31) 

where /(A) behaves as the zero in the addition and the 
identity in its multiplicative form. The problem of this 
definition is that G is in the corresponding (directed) 
cycle space but not necessarily in C{A). To avoid this 
unwanted situation, we propose the following projected 
addition for G in eq. (|30l): 



(32) 



Gi • G2 = G for G e £(A), 
Gi-G2 = G2 for G^£(A). 



In the second line, we call Gi is non-additive to G2, or 
vice versa. This projection drastically changes the al- 
gebraic nature of the operation. In fact, it breaks the 
commutativity: Gi • G2 7^ G2 ■ Gi, and the associativity: 
(Gi • G2) • G3 ^ Gi • (G2 • G3). So, the ordering of the 
operations matters in this case. It should be mentioned 
here that we can create some different definitions of the 
addition simply by different projections of G ^ 'C(A) onto 
£(A). As a starting point, we stick to the above defini- 
tion and will mention another possibility later. Note also 
that, precisely speaking, SG is either a candidate for a 
fiuctuation or a fluctuation which sometimes does noth- 
ing, because of the projection. 

If the configuration at the time Iq is given by G(to), 
the configuration G{to + t) at the time {to + t) can be 
expressed by a series of successive additions of 5G to 
G{to): 



G{h +t) = 



t/st-i 
Y[ SG{to + kSt) 



fe=0 



G{to). (33) 



The multiplication in the r.h.s. means the multiple ap- 
plication of the loop additions, OiLi = ■ ■ ■ A2 ■ Ai, 
and the addition should be done from right to left. The 
factor expresses nothing but a configurational and kinetic 
path between two configurations. 



Finally, we complete the definition of the stochastic 
fiow model by stochastically giving 5G in the same man- 
ner as in the previous section. Namely, assign the "po- 
tential" probabilities of a particle remaining at the same 
position and of a particle at i being kicked toward one 
of the nearest site j by p and x{i;j), respectively. For 
the isotropic distribution, x{i;j) can be simplified to 
x{i) = {1— p) / deg{i). Then, project the velocity changes 
onto the space of the oriented loop configurations C{h), 
and add it to the current configuration in the end. The 
degrees of freedom of 5vi and Vi are same as {deg{i) + 1), 
as long as an additional upper bound for 5vi is not im- 
posed. So, their spaces differ only by the meanings of p 
and x{i\j). For simplicity, we neglect this difference and 
distribute the probabilities isotropically. Note that the 
addition rule implicitly introduced the upper bound for 
5vi such that \5vi\ < 2. 

Now, we have formulated how the fiuctuations are 
given and added to existing flows. The propagator with 
the initial configuration Gi(io) can now be written down 
by: 



Zf"'^{x;G,,to;t + to) 



t/st 

SGt^steC(A) 5Gie£(A)J = l 

(34) 

where the probability P(G) of having the loop G in A is 
the same as in the stochastic Markovian model. Note the 
process of this model is also Markovian. 



IV. SIMULATION RESULTS 

In this section, we perform numerical calculations on 
the proposed models and show their statistical proper- 
ties for further investigation, some of which reveal very 
interesting features of the models. We restrict the lattice 
A to be as simple as the finite rectangular honeycomb 
lattice of roughly N x N hexagons for odd N, labelling 
it by A;i„ . Precisely speaking, we define the lattice A/i„ 
by piling {N — 1) and N aligned hexagons alternatively 
up to N rows so that the lattice becomes left-right sym- 
metric (Fig 13]). For odd N, the lattice becomes actually 
top-bottom symmetric as well, and therefore symmetric 
under discrete rotations. In the case of a finite lattice 
with boundary, there are the boundary sites and bonds 
which face the outer region of the lattice on its planar em- 
bedding. Among such sites, some have smaller degrees 
from those in the bulk. In the case of the honeycomb lat- 
tice, there are degree-two boundary sites. Accordingly, 
we can assign two different probabilities: x — {1 — p)/3 
in the bulk and Xb = (1 — p)/2 for the boundary sites 
of degree two. We choose the highest mobility p — 0, 
X = 1/3, and Xb = 1/2 for numerical calculations. 

As two models share the oriented non-intersecting loop 
configurations as an input to the system, the input data 
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FIG. 3: The illustration of the rectangular honeycomb lattice 
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has been generated in common for a million of time steps 
for various sizes, N = 11,21,31,41,51,61. Below is the 
basic information of A/i„ : 

dimViAh^) - (iV + l)(2iV + l)-4, 
dimVbiAh^) = 4iV, 

N+1 



2 2' 



F 



(35) 



for odd iV > 1. Vf, is the set of all the boundary sites of 
degree two. The number of faces is denoted by F. 



The statistics of the input data and the 
stochastic Markovian model 



We first plot the probability Pj(A/,„) of the no- loop 
configuration, or of the total frustration, for different 
sizes in Fig.Q. We have observed that it decays approx- 
imately exponentially with respect to the area of the lat- 
tice. Here, we propose the approximate functional form 
of Pf for general shape of the honeycomb lattice with F 
connected faces: 



(36) 



Note the form may not well approximate if A takes an 
extreme shape such as a lattice with small F or of aligned 
faces. Fitting the data in Fig|4]with the above function 
actually leads to no ~ 347.129 ± 10.53. We have also cal- 
culated Pf by hand in the case of the honeycomb lattice 
of 2 piles of 3 hexagons, A/j^ 3 , and confirmed that the 
extrapolation from A/i^ 3 gives rip ~ 355.089. This re- 
sult further indicates that the system gradually changes 
its nature between frustration dominant surface to loop 
dominant surface as the size changes around ^ 19. 

Next, we plot the number of events for the total length 
of the loops 1{G), which is equivalent to the number of the 
moving particles in one time step (FigjS]). It is remark- 
able that the leading contribution is not necessarily from 



FIG. 4: The probability of the no-loop configuration for vari- 
ous sizes in the log scale in the y-direction. The x-axis stands 
for the area of the lattice counted by the number of hexago- 
nal faces in A. There can be obviously seen a scalability. The 
error bars are given at the 95% confidence level. 



o 



c 
> 

LU 





life 










11x11 - 
21x21 - 
31x31 ■ 
41x41 
51x51 
61x61 ■ 












- 

□ 






,* 

\: ■- 














* Irs 

r 


i ■■ 
^ i 




■ 

as 


.■ 








■. % « 

1 ! iu ■;■ 

■ : lii II 






■■ If 

'm 

M 


Q 
















» : 
it 




il 

\ ^ 









20 40 



100 120 140 161 



Total length of the loops (in ht) 

FIG. 5: The distribution of the total length of the loops 
generated in one time step. The a;-axis is the length in the 
unit of the lattice spacing. The y shows the number of events 
in a million trials. Periodic oscillations seem to suggest the 
multiple entries of the same type of loops, e.g., two length-six 
loops at once contributes to the length-twelve event. 



length-six loops, and neither is the sub-leading contribu- 
tion from length-ten loops. For instance, on A/i^^, the 
leading and sub-leading contributions are from length-24 
and length-30 loops, resp. We have confirmed that the 
proposed form of P/ can explain nicely the A^-dependence 
of the probability of length-six loops. We also report the 
proportionality of the mean loop length I w.r.t. the area 
of the lattice: 



I - 0.01763F, 



(37) 



for F » 1. 

As for the loop number distribution, we find that they 
obey the Poisson distribution in an accurate way with 
the mean a approximately proportional to the number 



8 



of the faces F: 

a ~ 0.002856F. (38) 

The relation of the mean with the area of the lattice is 
not so trivial, suggesting that there is a relatively stable 
value of the loop number per unit surface as well as that 
of the total length of the loops per unit surface. Together 
with the average absolute value of the velocity and other 
basic numbers of A/i„ , we hst them in Table H] The very 



iV 


dimV 


dimE 


F 


/ 


a 


1^1 


11 


272 


386 


115 


2.01001 


0.326687 


7.38974 X 10'^ 


21 


942 


1371 


430 


7.558724 


1.226437 


8.02412 X 10"^ 


31 


2012 


2956 


945 


16.625344 


2.695635 


8.26309 X 10"^ 


41 


3482 


5141 


1660 


29.236702 


4.737545 


8.39653 X 10"^ 


51 


5352 


7926 


2575 


45.336774 


7.346373 


8.47100 X 10"^ 


61 


7622 


11311 


3690 


65.068168 


10.540758 


8.53689 X 10"^ 



TABLE I: The basic numbers of Ah„ for various N . 



low rate oi\v\ means that each element of the membrane 
behaves like a random walk with a high ratio of trapping: 
1 — |d|, and the system exhibits a very slow dynamics 
compared to the unit of time, 5t. 

In fact, we have confirmed that the model (I) shares the 
same property as the random walk in a cage in terms of 
the time-averaged mean square displacements (TAMSD): 

( xi{tf)TA = fdt' \xi{i' + t)- Xi{t')\\{Z<^) 

where we have used the discrete version of the above 
with ti = Q and = 10® in the unit of 6t. To show 
the difference from the independent random walks, we 
further obtain the time-average of the spacial two-point 
correlation functions, ( with i fixed at the centre 
of A^,, (FigUl): 

{\v.\\v,\)^-^ r dt'{\vM\vM)- (40) 

J-e Ej Jti 

As you see, our model has the notable two-point cor- 
relation functions which are clearly distinguishable from 
those of the simple random walks. 

B. The statistics of the stochastic flow model 

In the numerical calculation of the flow model, the 
initial configuration at t = is set to be the no-loop 
configuration /(A) for all sizes. The system gradually 
evolves as time goes, and finally reaches a quasi-critical 
state where all observable quantities become relatively 
stabilised. Their movements are schematically illustrated 
in Figl7]in the case of the periodic A/^. Although we 
have chosen the particular initial configuration, the data 
for various statistics is basically taken a certain amount 



Correlation (log) 




FIG. 6: An example of two-point correlation function of 
one of which is fixed at the centre of the lattice 
A;i2j . The z-axis is the number of events in a milUon trials in 
the log scale. 




FIG. 7: An illustration of the lattice Uquid on a honeycomb 
lattice. The sequences of green balls show the movements of 
the particles. The lattice is periodic to the horizontal direc- 
tion. 



of time after t = except for the probability of the fluc- 
tuations. We count how many times the fiuctuations are 
actually added to the system from the beginning to the 
end. Dividing it by the number of the input trials, we 
plot them in Figl8]as the probability. 

For the size iV > 41, the system becomes very stub- 
born about its stability against the "potential" fluctua- 
tions. Particularly for N > 51, much fewer fluctuations 
are added in contrast with the high probability of gener- 
ating the "potential" fiuctuations: Pi{Ahj^) ~ 1. There- 
fore, we will not look into these cases in detail, where the 
configurations are expected to be nearly frozen after it 
reaches a quasi-critical state. Note that there is no 'crit- 
ical state' such that no fluctuation can disturb, because 
there always exists at least one fluctuation trivially addi- 
tive to the present configuration G. It is the inverse 
that transforms it to /(A) by the addition. 

Instead of plotting the total length of the loops, which 
is well stabilised in a quasi-critical state, we plot the dis- 
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I, the average number of siniultaneous_ loops a, and the 
average absolute value of the velocity \v\ are all shown in 
Table im The above averages include the data from the 
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N: the size of the lattice 
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FIG. 8: The average probability of the fluctuations in the mil- 
lion steps. The probability suddenly drops down in between 
iV = 31 and = 41. 



tribution of the size of the individual loops in Figl5] At 
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FIG. 9: The distribution of the size of the individual loops. 
There is an almost uniformly scaling region in the middle. 

each time step, we classify the loops by their perimeter 
lengths, and accumulate them w.r.t. the length. So, for 
example, if there is only one oriented one-loop of length 
six existing during a thousand of time steps, the number 
of events counted for the length-six loops is a thousand. 
The distribution is drawn as such for various sizes, and 
we have found neither sharp peaks nor obvious cut-offs 
for the loop sizes in the plot. The most probable loop 
is the smallest loop of length six, and it seems to play 
a leading role of the dynamics of the flow model. There 
is no scaling relation in the entire region though, there 
is an almost uniformly scaling region in the middle part 
of the distribution. Moreover, the slope of the scaling 
seems to change suddenly between iV = 11 and = 21. 
Therefore, we infer that the rise of the loop probability 
Pi (A) for 5G may affect the evolution and the size of the 
individual loops. 

The average loop size, the average length of the loops 



N 


Ave. loop size 




a 


\v\ 


11 


17.015 


194.85 


11.452 


0.71638 


21 


19.259 


690.72 


35.865 


0.73325 


31 


20.120 


1484.8 


73.797 


0.73798 


41 


20.327 


2587.2 


127.28 


0.74303 



TABLE II: The basic data of the stochastic flow model on 
Ah„ for A = II,-- - ,41. 



beginning, and therefore we neglected the relaxation time 
tn from /(A) to quasi-critical states. The relaxation time 
Tjv which will be shown shortly is not very long compared 
to tall — 10^(5i. So, it would not spoil the accuracy of the 
statistics much. In this flow model, the distribution of 
the loop number behaves rather Gaussian than Poisson. 
In addition, I and a no longer obey simple proportional 
relations to the area of the lattice. The former has a 
rather correlation to the number of the sites dimV, as 
\v\ — l/dimV is almost independent of the size A^. The 
value \v\ gives the ratio of the moving particles among all, 
while 1 — |f I is of the trapped particles. Nearly constant 
value of |^;| might suggest that the value is intrinsically 
universal for all honeycomb lattices. The average loop 
number a is not so stabilised compared to the average 
loop length I. It should be remarked here that the pro- 
jection in our addition (I32p would certainly change the 
probability of the fluctuations, but only its characteris- 
tic time scale. So, we claim that the above values and 
qualitative arguments are characteristic in this type of 
models, regardless of the addition rule defined earlier. 

It is the dynamical aspect of the model that the addi- 
tion rule makes different. In this stochastic flow model, 
the particles are basically trapped in the sites alike in the 
stochastic Markovian model. But, once they are involved 
in the streamlines, they keep moving until it is stopped 
by a fluctuation. Coarse-graining of the pathway actu- 
ally make itself look like a Brownian motion, but it is 
actually much more complicated. Let us quantify it by 
the TAMSD of eq. gO]). We picked up a particle situated 
at the centre of the lattice and calculated the TAMSD 
(FiglTO|. where the lower limit of the time is set by the 
relaxation time ti — t^: tii_2i ~ 10^ (5i, r^i^ii ~ 10^ (Si, 
so that at i ~ Tat the particles are randomised enough 
on A from the initial configuration, tg is set to be the 
end of the simulation: the million-th step. As one can 
see from the left side of Fig[lUl the particle starts moving 
at a certain rate of \v\, and diffuses like liquid, slightly 
faster than the Brownian motion. At around t ^ 10, the 
exponent v in (a;(t)^) ~ f reaches the Brownian one, 
and later decreases to the anomalous diffusion of < 1 
in a strange way. Finally, it approaches its asymptotic 
value. By the times the particles approach the asymp- 
totic values, we can read off the relaxation times, tn. 
Their orders of magnitude are roughly the same we set 
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FIG. 10: The time-averaged mean square displacements 
(TAMSD) of the stochastic flow model (II) on A^^, for vari- 
ous sizes in the log-log format. The a;-axis is the time lapse t, 
while the y-axis is the square distance in the unit of the lattice 
spacing squared. The dashed lines are all exactly calculated 
asymptotic values for various A''. 



for the initial time ti for TAMSDs. It can be said that 
the system is generally in an anomalous diffusion process 
where v > lin the short range, and v < lin the distance. 
There can be found a strange and interesting plateau as 
well in the middle of the curve for each size of the lat- 
tice, but we have not yet obtained any clear explanation 
of the plateau. We will further mention this in the final 
section. 

There is a closest model to compare with our model 
(II), which is the non-reversal random walk. In that 
model, a random walker can bring a little memory for 
its past path so that it will not come back to the mem- 
orised positions. Usually, the walker memorises only the 
previous position in order not to come back at once. The 
model is sometimes known as a polymer chain model with 
excluded volume. In our flow model, a particle movement 
is always associated with a ring of movements and cannot 
reverse immediately. This minimum regulation matches 
with that of the non- reversal walker. In addition, because 
the non-reversal walker is not bounded by any loop-like 
relation to neighbours, we deduce that it gives the upper 
hmit for the TAMSD of our model. 

We define the stochastic flow model (II) and the non- 
reversal random walk on the lattice A/i^^ , and compare 
them in terms of their TAMSDs in FiglTTJ Notice that 
the constant trapping ratio is equipped with the non- 
reversal walkers to keep the pace of his first step with 
the flows in our model. 

As for the ergodicity, we claim that the ergodicity 
weakly holds in our model (II) of the lattice liquid. In 
order to illustrate this, we have calculated the asymp- 
totic value of the MSD with some assumption as follows. 
Assume that, in the limit of t — >■ +oo, the time corre- 
lation of every particle on A vanishes so that a particle 
situated at any of the sites on A at a certain time to is to 
be drifted to every site with equal footing in the infinite 
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FIG. 11: A comparison of TAMSD between different models 
with flows on Ah2i ■ The line for the non-reversal random walk 
is given with the trapping time consistent with the original 
model, and the asymptotic value for Ah2i is shown. 



future. Then, the asymptotic value of (a;/(t)^) can sim- 
ply be given by the average of the square distances of all 
the relevant pairs of the sites on A. On the rectangular 
honeycomb lattice A;j„ , we derived its exact formula and 



the leading terms are given by: 



>^o) = lim {xiity 

t^-\-oo 



1 



{7N^ + 1QN) +0{1). (41) 



In the case of iV = 11, 21, 31, 41, 51, the exact values are 
up to first five digits: 



116.75, 409.04, 876.45, 1518.9, 2336.4. (42) 



See Appendix ICl for the exact form of an arbitrary N. 

The time-averaged MSDs in FigfTUl [TT1 approach the 
calculated asymptotic values (jl^ and are stabilised 
around them as the time grows. Therefore, we can state 
that our models are surely ergodic in this sense. 



V. CONCLUDING REMARKS 

In this paper, we have proposed and formulated a fam- 
ily of models, called the 2D lattice liquid models, as one 
sort of the simplest models of single-component, single- 
layered, incompressible fiuid particle system on a mem- 
brane without explicit viscosity and density fluctuation. 
The model is formulated as a particle system on a lattice, 
but the configurations are found to be expressed by the 
sets of oriented non-intersecting loops, due to the restric- 
tion on the occupation number and the upper bound for 
the velocities. Following the formulation, we have con- 
structed two models of the lattice liquid. One of them 
is the stochastic Markovian model (I) which is found to 
be another manifestation of the vicious random walks in 
its densest phase and the DLL model in two dimensions. 
We have dealt with the model on the isotropic lattice. 
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and derived a few basic quantities of the model. We 
have also shown that the difference between our model 
and a collection of non-interacting simple random walks 
by distinct two-point functions. The other model is the 
stochastic flow model (II) where the set of rules are given 
in eqs. (123 [301 [32]). The first rule ^ is on the stability 
of the existing flow and is written in the Euler formu- 
lation. We have utilised the notion of the oriented loop 
configurations C{K) and formulated the model. In order 
to interpret the fluctuation as a direct consequence of the 
kicks by the environment, it is worth rewriting the rule 
in the original Lagrangian formulation. In this paper, we 
have focused on the construction of the model itself, and 
neglected this point. 

In numerical calculations, we have realised both mod- 
els on the rectangular honeycomb lattice Khn of various 
sizes in one million time steps. From the input data of 
the models on A/j^^ , we have found that the probability 
of the no- loop configuration, P/(A;i„), has a scaling re- 
lation with the area of the lattice, and we have proposed 
the functional form for a general size of the hon- 
eycomb lattice as an approximation. We briefiy report 
here that the probability of one-loop length-six configura- 
tions and that of length-ten can also be explained by the 
proposed form of Pf. This might suggest that the form 
of Pf with other additional criteria can explain how the 
distribution of the total length of the loops takes such a 
peculiar form in Fig[51 The area laws of the average loop 
length I and the average loop number d are also to be 
explained. We do not know the answers yet. In the nu- 
merical calculations of the flow model (II) , we have shown 
the probability of how many fluctuations of 5G had been 
added in a million trials. As N goes over 50, the system 
seems to be stuck against the projected addition. This 
reduction of the fluctuations is entirely due to the defl- 
nition of the projected addition. Because of the second 
line of the projection rule (15^ . if Gi is huge enough as a 
set, many elements of £(A) become non-additive. Even 
so, practically speaking in quasi-critical states, there are 
about 30% vacant faces. Therefore, the system cannot be 
entirely frozen, changing only the time scale for its dy- 
namics. It is the MSD on which the time scale matters. 
Especially, the emergence of the plateaus in FigfTUl might 
be because of this strong restriction of the fluctuations. 
We then expect that softening the projection in the ad- 
dition may change the intermediate plateau part of the 
MSDs, and this point is to be explored more elsewhere 

There are many related statistical models. Particu- 
larly, the vicious walk and the DLL model are the closest 
ones to the model (I). The vicious walk is the model 
of mutually avoiding multiple particle system so that its 
densest phase inevitably coincides with ours with the sim- 
plest stochastic process. The DLL model is the model of 
cooperative particle rearrangement for supercooled liq- 
uids and polymer melts, usually equipped with volume 
fluctuations and their temperature dependence. In both 
cases, the dynamics is governed by completely indepen- 



dent time steps, and they surely overlap with this sort of 
the models among the family of our models. Especially 
in the DLL model, the dynamics is usually equipped with 
the temperature dependence of a;(«), so it would be in- 
triguing to investigate another manifestation of the tem- 
perature dependence and multi-component systems (T2j 
in both contexts. In addition, the loop gas model is one 
of the closest models. Letting all Pf be one, one auto- 
matically obtains the partition function of the loop gas 
model from that of our stochastic Markovian model of 
the lattice liquid. The loop gas model can be described 
in the 0{N) spin model on the lattice by expressing all 
the loop graphs in terms of 0{N) spin variables [l3|. To 
avoid complications, we have not included this algebraic 
description of our models. Also, in the presence of the 
Pf less than unity, the loop probabilities of both models 
become different. We can quantify them by the use of 
number series known as the derangement number and its 
related series. To invoke the expression, the models must 
be first written in the form of algebras. We have put aside 
this procedure, and therefore, to be given elsewhere j7l| . 

We are afraid that the constructed models are not so 
applicable and not ready to be compared to experiments, 
except for a few cases in a very different context. Because 
we used the honeycomb lattice as an explicit example, the 
current model might be able to be applied directly to a 
particle system of the same situation: a bigger sort of 
particles form a most densely packed crystal sitting at 
the centres of the faces while smaller particles moving in 
between them. Because of the set-up, the MSD of the 
flow model actually resembles that of the super-cooled 
polymer melt flS^. Since their dimensions are different, 
we consider this happens by accident though. 

In order to apply 2D lattice liquid models to soft- 
material membranes directly, one should remove or re- 
lax some of the assumptions we used in this paper for 
convenience. Namely, the assumptions are such that the 
system must be single-component, single-layered, incom- 
pressible, non-explicitly-viscous, less density fluctuating. 
Also, we have not used interactions between particles. 
Note that they are not the requirements nor the defl- 
nitions of the 2D lattice liquid models. Among them, 
double-layer can easily be embedded just by folding a 
layer into two, but the interactions between two layers 
are not obvious. The simplest one is to introduce the 
interaction between layers only at the common bound- 
ary of two layers, possibly as a boundary condition. Our 
construction itself is ready to be extended to the multi- 
component case where there are more than one parti- 
cle species and the interaction between different species 
should be introduced as well. The simplest one will be 
the Ising-type nearest neighbour interaction with which 
one can observe the phase separation according to some 
temperature. This might also turn on the effective vis- 
cosity in the model. Note the concept of the temperature 
can naturally be installed by letting a;,p be a function of 
T, for example, by a; = e~^/'^ / deg and p = 1 — x ■ deg. 
Thermodynamic behaviour will depend on the manifes- 
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tation of the functional forms of x and p. Some of the 
above proposals will be studied elsewhere [ill 

Another very interesting yet possibly very hard model 
to solve is the fluid model with memory effect. As we 
mentioned in the manuscript, the "potential" probabili- 
ties can potentially carry the information from the past 
so that the past partially control the forthcoming event 
through the probabilities. Here, we have only dealt with 
Markov processes in both models, and the problem is 
beyond our scope of this manuscript. 
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Appendix A: On the space of the loop configurations 

To begin with, let us define the (oriented) loop config- 
urations on A. For a given lattice A, a one-loop graph, 
or configuration, is defined by a sequence of sites and 
bonds with no repeating sites and no ends. The ob- 
ject is the same as simple cycles in graph theory. An 
oriented one-loop graph is defined in the same way and 
is given by a one-loop graph with an orientation of the 
bonds representing its sequence. The space of the non- 
oriented non- intersecting loop configurations £(A) on A 
is the whole set of the non-intersecting combinations of 
such one-loop graphs. The space of the oriented non- 
intersecting loop configuration £(A) is those of the ori- 
ented one-loop graphs. These spaces can be decomposed 
into (oriented) fc-loop subspaces, or can be defined by 
them. 

Let Ci be a one-loop graph on A for i G (1, • • • , k). 
Then, a fc-loop graph Gk is defined by: 

Gk = Uti A, s.t. A n = for i ^ J. (Al) 

Collecting all such fc-loop graphs in £(A) constitutes the 
fc-loop subspace £fe(A) of C{A): 

£fe(A) = {Gfc}, £fe(A)c/:(A). (A2) 

When A is finite, there exists the maximum number of 
co-existing non-intersecting loops on A. Denoting it by 
kmax, one can write down the decomposition as: 

/:(A) = ©^— /:fc(A), 

s.t. £,(A) n £j(A) =0 for i 7^ j. (A3) 
In the same manner, one can decompose the oriented 



space C{A) into the oriented fc-loop subspaces: 

£k{A) = [Gk}, A.(A)C£(A), 
£(A)-e^:— 4(A), 

s.t. A (A) n Cj{A) =0 for i 7^ j, (A4) 

where Gk is an oriented fc-loop graph in C{A). If we intro- 
duce the function D, which counts the number of elements 
of a set, then the uniqueness of the no- loop configura- 
tion in both spaces is trivially given by Co{A) = £o(A) 
and n{£o{A)) — 1. We label this unique configura- 
tion by /(A). Similarly, n{Ci{A)) = 2n{Ci{A)) due to 
the orientations. One can further observe the relation: 
r!(£fc(A)) = 2'= niCkiA)) for < fc < fc^^x- 

The algebraic structure of the space C{A) directly cor- 
responds to the physical relations between the states on 
A. However, even for a very basic question of how many 
configurations are in the space C{A), we have not known 
any clear answer yet for an arbitrary size of the lattice. In 
fact, it is shown that counting the number of self-avoiding 
walks (SAWs) for fixed ends on a given graph is #P- 
complete [l3| . #P is a class of the problems for counting 
and, roughly speaking, they are expected to be as hard as 
A^P-complete problems such as the prime number decom- 
position of an arbitrary large number. Therefore, one can 
expect that counting the number of the elements of C{A) 
is as hard as A^P-complete problems as well, because a 
non-intersecting loop is nothing but a self-avoiding walk 
ending at its start |15| . 

Below, we illustrate the structure of the spaces in terms 
of the cycle space in graph theory. We restrict ourselves 
to the planar cases unless otherwise stated. When A is 
planar, i.e., A can be embedded on a plane with no bonds 
crossing, the space C{A) can be spanned by the set of the 
smallest loops which are the perimeters of the interior 
faces of A. In graph theory, there is the so-called cycle 
space cycle{A) which is generated by all the simple cycles 
and whose basis is given by the fundamental cycles. The 
addition generates the cycle space from the basis, and 
makes it the vector space of Zf-^+\ where E and V 
are the number of the bonds and the number of the sites, 
respectively. The concepts of C{A) and the cycle space 
cycle{A) are very close to each other, but not always the 
same because of the presence of the intersecting loops in 
the cycle space. In fact, there is a necessary and sufficient 
condition for their equivalence that the degrees of the 
sites must be less than four everywhere: deg{i) < 4 for 
all i e C{A). Otherwise, the addition allows intersecting 
loops. Therefore, 

C{A) C cyde(A), (A5) 

in general, and C[A) is not necessarily a vector space but 
can be spanned by the set of the smallest loops. 

When A is a planar honeycomb lattice of a finite size, 

£(A) = cycle{A) ~ Zf-^+\ (A6) 
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and 

0(/:(A)) = 2^-^+1 = 2^, (A7) 

where F denotes the number of the interior faces and the 
second equahty is given by Euler's formula on a planar 
graph. Its basis can be given by the set of the smallest 
loops. Note that, in general, 

n{C{A)) < 2^-^+\ (A8) 

In the case of the oriented loop configurations C{A), the 
relation is not as simple as above. We can only provide 
an estimate for the number of the elements to close. In 
the case of kmax — ~^ ? 

2^+1 - 1< r!(£(A,„,3„J) < 

2l^-2^(i^ + l) + 2i^ + l. (A9) 

A more refined relation will be given elsewhere [llj . 

Appendix B: The definition and the basic properties 
of the probability functions Pf and Pi 

In this section, we define the probability functions of 
the no-loop configuration and of having (oriented) loops 
on the given graph, Pf and Pi, respectively. That is to 
say, the complementary functions satisfying: Pf + Pi — 1. 
Upon their definitions, we will show some of the basic 
properties of the probability functions. Note that the 
word 'loop graph' indicates the same object as the loop 
configuration in the original context, but the latter em- 
phasises an outcome of some stochastic process. 

On the loop graphs as events of the stochastic process, 
we define the following probability functions. Let P{Gk) 
be the probability of having the loop graph Gk in A. Let 
P/(G) be the probability of having any loop graph on 
the graph G A. Let Pf{G) be the probability of the 
graph ff^ G A being completely frustrated, i.e., having 
no loops. Here, the graph is an extended graph of G 
such that it additionally includes all the bonds adjacent 
to the sites in G, in order to precisely explain and esti- 
mate the stochastic nature of our model. Introducing the 
complementary graph £ = A — £ to a loop graph £, the 
probability function P{C) of having the loop graph C in 
A is given by the product of the probability Pi (£) of the 
loop configuration C on the graph C, and the probability 
Pf (£) of all the other sites being frustrated: 

P{C) = PiiC) ■ PfiC). (Bl) 

By definition, the following equalities and inequalities 
hold: 

Pi{A) = Yl ^(G*^-), 

Gfc>oe£(A) 

PfiG) - 1 - Pi{G), Pf{%) = 1 - Pi{%) = 1, 

< < 1, 0<P/(G)<1. (B2) 



In the third line, for finite G and A, the equalities 
P;(G) = 1 and P/(G) = only hold when G or A is lin- 
ear or the assignment of the "potential" probabilities of 
the particles' moving are unnaturally given in the non- 
stochastic limit. Therefore, we assume < Pi{G) < 1 
and < Pf{G) < 1 for finite size G and A. 

Letting C be the naive projection of the oriented loop 
graph C by removing the orientations, and assuming that 
the different oriented loop graphs which share the same 
projected loop graph C are equiprobable, then the map- 
ping of the oriented loop graphs by the probability func- 
tions become apparent as: 

P(Gfc) = Pi{Gk) ■ Pf{Gk), 

P{Gk) = 2^-p(Gfe), 

Pi{Gk) = 2''P{Gk), 

PfiGk) = P/(Gfe). (B3) 

Gfc is not a subgraph of A so that the complementary 

graph is defined by the projected graph Gk- Gfc = A — 
Gfc = Gfc. One can interpret that the same projection is 
done for the Pf as well. 

Because the frustrations in disconnected graphs are in- 
dependent, the following relation holds: for two discon- 
nected graphs, Gi,G2 £ A: 

PfiGi n G2) = P/(Gi) • P/(G2) for Gi n G2 - 0.(B4) 

Note that, using the above definitions and successive ap- 
plication of the relations, one can obtain P(Gfc) in a 
nested expression of the loop probabilities only. 

Appendix C: Asymptotic value of the mean square 
displacements (MSD) 

1. Set-up and derivation of the exact asymptotic 
values of MSD 

In this section, we calculate the asymptotic value of the 
mean square displacements, ( R^{t)), which is the average 
of ( x{tY) over /, in the limit oft ^ 00. Before going into 
the actual calculation, we fix the notations and introduce 
some functions to facilitate it. 

Number the position on A starting from to M— 1, and 
assign a pair of integers tothep-thsite, p = {px ,Py). The 
position of the p-th site, (xp, ijp), on the two-dimensional 
plane can be given by, for the N x N honeycomb 
lattice: 

Px = p mod col, py — [pjcol], 

xp^^Px, y,^\py + \{-iY^^'^^\ (CI) 

where M = dim{V{A)) = {N + 1){2N + 1) = 2N^ + 
37V -1-1, col = 2A^-|-1, and [A\ denotes the greatest integer 
which does not exceed A. px takes the value from to 
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{col — 1) = 2N, while Py = 0, . . . , N . p^^y are all integers. 
The lattice spacing is taken to be a unit length. The 
square distance between p-th and q-th positions are: 



more than one bond. For instance, regardless of the value 
of N, the following two sites should always be excluded: 



q,p-q 
3 



= {xp - Xqf + {yp - ygY 



(7 -{(0,0), (27V, 0)}. 



(C5) 



{{Px - Qx f + HPy - Qy)'^} 

+ l{l_6(-l)f==+P«(p,-g,)}(l-(-l) 



(Px-9i) + (Pa 



Note that we used the following identity: 1 = (—1)^*^ for 
Vfc G Z. Replacing p — q by (fc, I), one can rewrite this 
as R, 



q.Xk.l)- 



(C3) 



Besides, when N is odd, the following two sites are to be 
taken into account in the same manner: 

q^mN),{2N,N)}. (C6) 

(C2)rhey are connected only one nearest neighbouring site, 
so that the particle cannot reach these four points via our 
flow dynamics. Therefore, the square distances related to 
these points do not contribute to the above calculation. 
Substituting (jCSp into (IC4|) and subtracting such unnec- 
essary terms from the sum, we obtain for odd N: 

(Rl) 



where C(q, /) = (! + 6/(— 1)''^+'''') . From this expression, 
one can verify that i?^ ^j, is integer. 

Assume that, in the limit of t — >■ +oo, all the pairs 
of the sites are independent to each other so that the 
asymptotic value of (i?^(i)) can simply be given by the 
average of the square distances of all the relevant pairs 
in M X M combinations. Namely, the value leads to the 
expression: 



{{N + 1){2N+1) - 4 + (1 + 
-{N + lf{2N + lf{7N^ + ION + 1) 



+ -(iV + l)^(2A^+l) 



Ri 



\im{R'{t)) 

t—¥oo 



-2 {UN^ + 29N^ - + N + 1) 
-i{6(iV+lH2A^+l) + l}(^i±i^ 



(M2),, 



levant 



N col-1 N-qy col-l-q^ 

E E E E 



relevant 

(C4) 



;^14A^'^ + 29iV^ - 16N^ -2N + 2) 



;{7N^ + ION), 



l + (-l) 



N 



(C7) 



Here, the word "relevant" means that some sites at the up to 0(1). In the cases of = 11,21,31,41,51, the 
boundary cannot contribute to the above, because our exact values are given in the manuscript up to first five 
fiow dynamics can only involve those points which have digits. 
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